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Abstract 

The Passive Aerodynamically Stabilized Magnetically Damped Satellite (PAMS) will be deployed from the Space 
Shuttle and used as a target for a Shuttle-mounted laser. It will be a cylindrical satellite with several comer cube 
reflectors on the ends. The center of mass of the cylinder will be near one end, and aerodynamic torques will tend 
to align the axis of the cylinder with the spacecraft velocity vector. Magnetic hysteresis rods will be used to 
provide passive despin and oscillation-damping torques on the cylinder. 

The behavior of the hysteresis rods depends critically on the "B/H" curves for the combination of materials and 
rod length-to-diameter ratio ("1-over-d"). These curves are qualitatively described in most Physics textbooks in 
terms of major and minor "hysteresis loops". 

Mathematical modeling of the functional relationship between B and H is very difficult. In this paper, the physics 
involved is not addressed, but an algorithm is developed which provides a close approximation to empirically- 
determined data with a few simple equations suitable for use in computer simulations. 
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II. Motivation For Analysis 


A few years ago, a proposed Gravity and Magnetic Earth Surveyor (GAMES) spacecraft design included a small, 
totally passive subsatellite which would be released to fly at controlled distances behind the main spacecraft. It 
was to carry a single laser-reflecting comer cube which ideally would be pointed directly at GAMES and provide 
a target for a laser-ranging experiment. 

A novel subsatellite passive attitude stabilization method, using a combination of aerodynamic and magnetic 
torques, was proposed by Dave Skillman and Jim Abshire of the Goddard Space Flight Center (GSFC). The work 
reported here concerns the behavior of magnetic hysteresis rods which were to be used for subsatellite despin and 
oscillation control. 

Although a GAMES new start did not materialize, the Project recommended, and received approval for, a Shuttle 
test flight of the aero/ mag stabilization concept. Preparations are currently underway for the flight of the Passive 
Aerodynamically Stabilized Magnetically Damped Satellite (PAMS) in 1996. 

III. Development of the Boundary Curves 

Magnetic hysteresis rods (long, slender, cylindrical samples of a permeable material) become magnetized when 
exposed to an ambient magnetic field. The "magnetizing field" (H) is the component of the external field parallel 
to the axis of the cylinder, and the "magnetic flux density" (B) developed in the material produces a magnetic 
moment parallel to the axis. The interaction of the magnetic moment and the external field produces a torque 
which is sometimes used to despin satellites. 

For the type of materials used, there is a significant phase lag between B and H and most physics books show 
plots of "B/H curves" or "hysteresis loops". The function B(H) is extremely nonlinear and multi-valued. 

There is no physics addressed here. The mathematical model suggested merely attempts to reproduce empirical 
B vs. H behavior in a form suitable for use in computer simulations of dynamic systems which employ magnetic 
hysteresis rods. 

We begin with the familiar function y = tan(x) shown on the left below in Figure la. Interchanging the axes 
produces the not-as-familiar function y = arctan(x) as shown below in Figure lb. 
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We now let x = kH where k is a selectable constant 
and H is an independent variable. We also let y be 
proportional to a dependent variable (B) and wish to 
have the values of B limited between "saturation 
levels" -B s and +B S . The appropriate relationship is 
given by 


In terms of the new variables, we now have the 
function 

B = — B s tan -1 (kH). (2) 

71 

Figure lc shows this function, where k has been 
arbitrarily set to 5.92 and B s to 9872. 



Figure lc B = — B s tan ^kH) 
71 


Figure Id shows the result of shifting the previous 
curve to the left by an amount H c arbitrarily set to 
0.135 units. Note that, where this curve crosses the 
ordinate axis, we let B = B r . The equation for this 
curve is 

B = — B s tan -1 [k(H + H c ) 1 . (3) 

71 1 J 

For H = 0 then, we have 


B r =-~B s tan" 1 [k(H c )]. 
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Figure Id B = — B s tan 1 [k(H + H c )] 
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Figure le shows the function 

B = |B s tan- 1 [k(H-H c )], (6) 

obtained by shifting the original curve to the right 
by the same amount H c 


Finally, in this construction sequence, we 
superimpose the last two figures to obtain Figure If 
which bears a close resemblance to the "magnetic 
hysteresis loops" shown in the Electromagnetics 
chapter in most elementary Physics books. Usually, 
H is called the Magnetic Field Strength and B is 
called the Magnetic Induction. The units involved 
are discussed in an Appendix to this report. Note 
that the left boundary curve crosses the ordinate axis 
at B=B r . In textbook B/H curves, this value is called 
the Remanence. The right boundary curve crosses 
the abscissa at H=H C . This value is called the 
Coercive Force. 



Figure le B = — B s tan ^(H-Hj] 



IV. Initial Magnetization Phase and Behavior Inside the Boundaries 

Most Physics texts mention an initial magnetization phase in which B(0)=H(0)=0 and H is gradually increased in a 
positive direction. Figure 2 shows a typical "S"-shaped curve showing B first slowly increasing, then rising more 
sharply, then asymptotically approaching some upper limit. 

Here, we assume that the slope dB/dH of this curve and that of any point within the boundary constraints 
depends on the horizontal distance between the current point (H,B) and the boundary curve which the point is 
moving away from, in this case, the left side boundary. As an example, suppose that the slope is zero coming off 
the left boundary and approaches the boundary slope as it approaches the right boundary. Note from the 
previous section that the two boundary curves are always 2H C apart in the horizontal direction and define the 
fraction 
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f _ H-H L 

2H C 


(7) 


where H L is the value on the left boundary curve 
corresponding to the current value of B. Since the 
left side boundary is given by 


then. 


B = -B s tan- 1 [k(H L+ H c )], 


tan 


H, 


^7C B ^ 

v2B s y 


-H c . 


(8) 


(9) 



Now, Let B' be the boundary curve slope (either boundary) corresponding to the current value of B. 
Differentiating the previous equation with respect to H, we find 

n 2 f n B ^ 

sec 

2kB c 2 


0 : 


B' 


's J 


or 


(10) 


B' = — kB, cos" 

71 
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Now, a simple function providing the desired behavior could be 

dB 


Another possibility would be 


dH 

dB 

dH 


= f B' . 


= fP B' 


where p is an exponent of the fractional distance f. 


( 11 ) 


( 12 ) 

(13) 


Actually, in the specific example addressed below, the slopes at the left boundary seen in the test data were small 
but non-zero. To account for this possibilty, a more general expression (with one more selectable constant) was 
used in fitting the data. We let 


dB 

dH 


[qo+(i-q 0 )f p ]B' 


(14) 


where q Q B' is the slope dB/dH near the boundary. 
V. Specific Application 


Some empirical B/H curve data was generously provided to the author by Frederick Mobley of the Johns 
Hopkins University Applied Physics Laboratory (APL). Tests were run on a long, slender rod made from AEM 
4750 material, a nickel-iron alloy. The magnetizing field was oscillated between nine different peak values 
ranging from +/- 0.04 to +/- 0.50 units. At the higher levels, boundary curve limits were clearly achieved. 

Examination of the boundary curves showed immediately that H c =0.135 units and B r =4240 units were reasonable 
approximations. Placing these values and H=0 in the left-hand boundary curve 



tan _1 [k(H + H c )] 


(15) 
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( 16 ) 


produces an equation showing a required relationship between k and B s/ i.e. 

4240 = | B s tan _1 [0.135k]. 

Another point on the left hand boundary was located at (-0.5,-7150). Using these values for H and B, respectively, 
produces a second requirement for k and B^, i.e. 

-7150 = — B s tan _1 [-0.365k|. (17) 

71 

These equations were solved numerically to find B s =9872. With this known, we can compute k = 5.925. The 
resulting boundary curve then becomes 

B = ^ (9872) tan -1 [5. 925(H ±0.135)]. (18) 


Figure 3a shows a reproduction of some of the 
hysteresis rod B/H test data received from APL. 
The outermost contour was obtained by cycling the 
magnetizing field (H) between the limits of ±0.5 
units and plotting the resultant steady-state B/H 
contour curve. The peak values of B observed were 
±7150 units. Points on the left-hand boundary curve 
were used to determine appropriate values for k and 
B s for the math model described above. 



The outermost contour of Figure 3b shows the math 
model output for the same ±0.5 units input. Note 
that the test data slope appears to have a fairly large 
linear range compared to that produced by the 
model. The arctangent model shows more curvature 
and flattens out a little more at the extreme values 
of H. 



Figure 3b Model Data 
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Figure 3c shows test data for a smaller range of H 
values, i.e. ±0.1 units. The outermost curve here is 
the same as the innermost curve of Figure 3a. This 
curve was used in the selection of the other two 
model parameters, p and q 0 Values of p = 4.75 and 
q 0 = 0.085 were obtained by a trial-and-error iterative 
procedure which matched a peak value of B of 637 
units and an ordinate crossing value of 275 units as 
observed in the test data. Figure 3d shows that a 
very close match was obtained. 


Figure 3c Test Data 



The four parameters selected above were then used 
in the model andB(H) contours were generated for 
the other ranges of H for which test data was 
available, i.e. ±0.4, ±0.3, ±0.2, ±0.15, ±0.08, ±0.06, and 
±0.04 units. All test data and model outputs are 
shown in Figures 3a through 3d. 
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Figure 3d Model Data 





VI. An Algorithm for Determining dB/dt 

Given B, H and dH/ dt, with k, and p predetermined constants, 

1. Compute H L/ the value of H on the left boundary curve corresponding to B. Since this curve is given by 
B = - B s tan" 1 [k( H L + H c )] , we find 

HL=TAN(PI*B/2/BS)/K-HC 


2. Compute B’, the boundary curve slope corresponding to the value of B. It may be shown that 


3. Find f. 


BP=2*K*BS/PI*COS(PI*B/2/BS)**2 

F=(H-HL)/2/HC 


If dH/dt is negative, measure f from the right hand boundary. 


4. Find q, 

5. FinddB/dH, 

6. Find dB/dt, 


IF (DHDT.LT.O) F=l-F 
Q=Q0+(1-Q0)*F**P 


DBDH=Q*BP 
DBDT =DBDH*DHDT 


VII. Appendix on Magnetic Units 

The use of units for the variables B and H has been avoided in the body of this report. The history of magnetic 
units is long and sometimes contentious and it was feared that including a discussion of this subject would 
introduce an unnecessary distraction in the presentation of the mathematical modelling of magnetic hysteresis. 

As mentioned in the body of the report, empirical B/H curves were received from Mr. Fred Mobley of The Johns 
Hopkins University Applied Physics Laboratory (APL) of Silver Spring, MD. The data referred to here was 
produced during testing of a rod made from a material known as AEM 4750 (a nickel-iron alloy) with a length of 
of 58 inches (1.47 meters) and a diameter of 0.1084 inches (0.00275 meters). The independent variable, the 
magnetizing field H, was expressed in oersteds and the dependent variable, the flux density B was given in 
gauss. 

In modem times, people have become comfortable with magnetic fields in tesla, magnetic moments in ampere- 
meters-squared and torques in newton-meters. A magnetic hysteresis rod in space will have a certain magnetic 
moment M and will interact with the Earth’s magnetic field B to generate a torque T, with the pertinent equation 
being the familiar cross product law T = M x B. 

In this Appendix, we attempt to bridge the gap between raw data in the form of B/H curves, with B in gauss and 
H in oersteds, and computer simulations of the performance of magnetic hysteresis rods in space. 

An immediate possible source of confusion is the fact that the letter B is now commonly used for the external 
magnetic field, where historically, the letter H was used. In early-space-age literature 1 and textbooks 2 , when 
magnetic materials were exposed to a magnetic field strength or magnetizing field H, a magnetic flux density or 
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magnetic induction B was developed in the material. The unit for H was oersted and the unit for B was gauss. 
These unit names honored two of the early researchers in the field, Hans Christian Oersted (1777-1851) and Carl 
Friedrich Gauss (1777-1855). 


With a magnetic hysteresis rod, the magnetic moment or dipole moment depended on the flux density and the 
volume of material with the governing relationship being m = BV / 4ji. Here the magnetic moment m is measured 
in pole-centimeters and the volume V is in cubic centimeters. Pole-cm is a unit probably familiar to most readers 
and many have probably used the conversion factor of 1000 pole-cm = 1 ampere-meter-squared. Why this is so is 
probably not very well known. 


The ampere-meter-squared unit is probably the simplest and most obvious magnetic unit. With a current in 
amperes flowing in a loop enclosing an area in square-meters, the magnetic moment is simply the current times 
the area. Now, the magnetic field in the vicinity of such a loop is of interest. The physics involved is expressed 
bv the Biot and Savart equation 

dH = I ds x r / r 3 . (19) 


In this differential equation, H is the magnetic field, I is the current, ds is an incremental length of coil and r is a 
vector from ds to the point of interest. Unitwise, with r and ds in centimeters, a current in abamperes produces a 
magnetic field in oersted. There are 10 amperes in one abampere. 


Figure A-l shows a circular coil of radius "a" 
carrying a current "i" (in amperes) and a point R in 
the plane of the coil at which the magnetic field is to 
be determined. For points far away from the coil, 
the following analysis shows that the magnetic field 
is given approximately by 


H = na 2 i/10 R' 3 
oersted "into the paper". 


( 20 ) 



The elemental length ds is in the 1-2 plane, with the components ds } - - a sin(0) d0 

and ds 2 = a cos(0) d0. 


The vector r (also in the 1-2 plane) has the components 


r t = R - a cos(0) 
and r 2 = -asin(0). 


( 21 ) 

( 22 ) 

(23) 

(24) 


The vector ds x r is obviously in the 3 direction with a magnitude of a[a-R cos(0)]d0. The quantity r 2 is given by 
R 2 +a 2 -2aRcos(0). The strength of the field is then obtained by integrating (from 0 = 0 to 2 n) the differential 


equation 


dH = i/10 a [ a - R cos(0) ] d0 / [ R 2 + a 2 - 2 a R cos(0) ] 3/2 . (25) 


Far from the coil, i.e., for R»a, the denominator in this equation may be approximated by R 3 [l-3(a/R)cos(0)] in 
the denominator or by R' 3 [l+3(a/R)cos(0)] in the numerator. This makes the approximate differential equation 

dH = a i/10 R- 3 [ a - R cos(0) ] [ 1 + 3 a/R cos(0) ] d0. (26) 
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(27) 


Of the four terms here, only two will contribute to the integral, so we can let 

dH = a i/10 R' 3 [ a - 3 a cos 2 (0) ] d0. 


Integrating this equation produces 


H = a i/ 10 R' 3 ( 2 ji a - 371 a ) = - na 2 i/10 R 3 . (28) 

Now, fictitious "magnetic poles" were commonly employed in the teaching of magnetics years ago. They came in 
two varities, positive poles and negative poles. The force between two poles was inversely proportional to the 
square of the distance between them and was measured in dynes when the distance was in centimeters. A dyne 
was the force required to accelerate a mass of one gram to a level of one cm/sec 2 . Forces were either repulsive or 
attractive, depending on whether the signs were the same or different. If a positive pole and a negative pole were 
placed one centimeter apart, they formed a unit dipole, with a magnetic moment of one pole-cm. Collections of 
p positive poles and p negative poles a distance "c" apart formed a dipole with a magnetic moment of m=pc. 
Now, a unit positive pole in the vicinity of this dipole would experience an attractive force toward the negative 
end and a repulsive force away from the positive end. The magnetic field (in oersteds) due to the dipole at this 
point was defined as equal in magnitude and direction to the net force vector (in dynes) acting on the unit pole. 


Figure A-2 shows the dipole described above and a 
point R at which we would like to determine the 
magnetic field strength. It will be shown below that, 
for R»c, and with R perpendicular to the dipole 
axis, the field is approximately H = pc/R 3 . Assume 
now that the fields due to the coil described above 
and the dipole described here are equal. We then 
would have 7ta 2 i/10=pc. Now, let A be the area of 
the coil in square meters. Then a 2 =10 4 A, and we 
have 1000Ai=pc. Therefore, one ampere-meter- 
squared = 1000 pole-cm as stated above. 



Figure A-2 Magnetic Dipole 

Initially, let R=Rr be a general point at a distance R from the center of the dipole, with r a unit vector. Let m be a 
unit vector from the positive pole to the negative pole of the dipole. Let P be a vector from the test point to the 
negative pole of the dipole and let Q be a vector from the positive pole to the test point as shown. We can write 
R+P=mc/2 and Q=mc/2+R. The force (in dynes) acting on the test pole will be pQ/Q 3 + pP/P 3 , where P and Q 
are the magnitudes of the vectors P and Q. 

Now, P 2 =PdotP=(mc/2-R)dot(mc/2-R) or approximately, with R»c, P 2 = R 2 - c mdotR. Obviously, we have 
introduced the operator "dot" to indicate the scalar product between two vectors. Similarly Q 2 is approximately 
equal to R 2 + c mdotR. The force then becomes, approximately, 

f = pQ [ R 2 + c mdotR ]" 3 ^ 2 + pP [ R 2 - c mdotR ]- 3 ' 2 (29) 

“ (P/R 3 ) [ Q ( 1 + 3/2 c/R mdotr) + P ( 1 - 3/2 c/R mdotr) ] 

« (p/R 3 ) [ P+Q - 3/2 c/R (P-Q) mdotr] 

= (p/R 3 ) [ c m - 3 c r mdotr ] 

® (pc/R 3 ) [ m - 3 mdotr r ] 
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Recalling that the force in dynes is numerically equal to the magnetic field in oersteds, this equation represents a 
concise expression for a dipole magnetic field. It is similar to that shown in Wertz 3 . 

At a point on the dipole equator (which is what we were initially interested in) the field reduces to simply pc/R 3 
as was stated above. 

Getting back to the computer simulation of magnetic hysteresis rod behavior, we have the following rule: 

Given a magnetic flux density (B^ in gauss . and a volume of material (V) in cubic cent imeters, the magnetic 
moment in modem (SI) units will be M = BY /4000 k ampere-meters-squared.. 

Another perhaps familiar expression from the magnetics section of physics books is B = pH, where p is the 
permeability of the medium. This equation converts H in oersteds to B in gauss. With the now obsolete (but still 
widely used) system of units, the permeability of free space was unity. Magnetic fields in space were then 
expressed in gauss and were numerically equal to magnetic fields in oersteds. This latter unit soon became rarely 
used and the Earth's magnetic field strength was stated in gauss (or milligauss) for some time. The unit gamma 
was also sometimes used in the early days of orbiting spacecraft with one gauss being equivalent to 10 5 gamma. 
When the emphasis on SI units arrived on the scene, "gauss" gave way to webers per square meter with the 
conversion factor of 10 4 gauss per Wb/m^. This unit name was selected to honor Wilhelm Eduard Weber (1804- 
1891). Some time later, this latter unit was renamed the tesla in honor of Nikola Tesla (1856-1943). This made one 
gamma equal to one nanotesla. 

A second rule to be followed then in dealing with B/H data in gauss and oersteds is: 

Working with magnetic hysteresis rods in space, the "H in oersteds" req uired for the math model would be 
obtained by multiplying the magnetic field in tesla by 10000. 
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